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ABSTRACT 

We discuss the intergalactic propagation of ultra-high-energy cosmic rays (UHECRs) 
with energies E > 10 18 ' 5 eV. We consider the propagation of UHECRs under the in- 
fluence of the energy-dependent deflection by a weak random magnetic field in the 
intergalactic medium and energy losses by photo-pion and pair production. We cal- 
culate arrival spectra taking full account of the kinematics of photo-pion production 
and the Poisson statistics of the photo-pion interaction rate. 

We give estimates for the deflection of UHECRs from the line of sight to the 
source, time delays with respect to photons from the same source, arrival spectra 
and source statistics. These estimates are confirmed by numerical simulations of the 
propagation in energy evolution of UHECRs. These simulations demonstrate that the 
often-used continuous approximation in the treatment of energy losses due to photo- 
pion production on the cosmic microwave background (CMWB) cannot be justified 
for UHECRs. 

We discuss the implications of these results for the observed flux of particles 
above the Greisen-Zatsepin-Kuz'min cut-off in the two of the scenarios that have been 
proposed for the production of these particles: continuous production in the large 
shock waves associated with powerful radio galaxies, or possibly large-scale structure 
formation, and the impulsive production at relativistic blast waves associated with 
cosmological gamma-ray bursts. 
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1 INTRODUCTION 

Recent observations of ultra-high energy cosmic rays (UHE- 
CRs) above 10 18 ' 5 eV (Bird et al. , 1994, Yoshida et al. , 
1994, Takeda et al. , 1998) suggest that these particles are 
of extragalactic origin. The indications are threefold: 



(i) The energy spectrum hardens, from J{E) oc E 3 - 2±01 
below 10 18 ' 5 eV to J(E) oc £r 2 ' 8±0 ' 3 above 10 18 ' 5 eV; 

(ii) Arrival directions do not seem to be clustered signif- 
icantly along the Galactic plane, although there is an indi- 
cation that they may be clustered along the supergalactic 
plane for E > 10 19 5 eV (Hayashida et al. , 1996; Stanev et 
al. , 1995). 

(iii) There is a suggestion that the composition changes 
from one dominated by heavy nuclei (iron group) below 
10 18 ' 5 eV, as expected for a population of Galactic origin, to 



one dominated by light nuclei or even protons above 10 
eV; 

At these energies the deflection of UHECRs in the ~ 
pG Galactic magnetic fields is too small to randomise flight 
directions within the Galaxy. In terms of the typical gyration 
radius of a particle with energy E and charge q — Ze and 
pitch angle i, pg = r g (_B)sini with r s (E) = E/ZeB, the 
deflection angle A8 inside the Galactic disk of size d ~ 10 
kpc is of order 



AO: 



d 



6° (ZB„ G /E 2 



(1) 



with E20 = -E/(10 20 eV). In practice, it will be much less 
due to the inclination of the particle orbit with respect to 
the disk plane and the fact that the Galactic magnetic field 
is strongly turbulent on kpc scales (for a review: Zweibel 
& Heiles, 1997), and will not act coherently over the whole 
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disk. Therefore, an origin of UHECRs within the Galactic 
disk should lead to a significant clustering of arrival direc- 
tions around the Galactic plane, unless the sources are in an 
extensive halo around the Galaxy. 

From a theoretical point of view, the fundamental limit 
on the maximum energy attainable in any electrodynamic 
acceleration mechanism relying on bulk plasma motions (i.e. 
waves or shocks) is (Norman, Melrose & Achterberg, 1995, 
Gallant & Achterberg, 1999) 

^max Zer s f3 s BR s . (2) 

Here j3 e = V/c is the wave or bulk velocity in units of the 
speed of light, F s = (1 — /3s)~^ 2 is the associated Lorentz 
factor, B is the magnetic field strength and R B is the linear 
size of the source. This limit strongly constrains any model 
for the production of UHECRs (Norman, Melrose & Achter- 
berg, 1995). Galactic sources, with the possible exception of 
radio pulsars, seem to be excluded. The value of -E max for 
these source falls well below the UHECR energy range. 

For the above reasons, most theoretical discussions of 
the origin of these particles have centered on extragalactic 
sources. Leaving aside for the moment the more exotic mod- 
els involving some fundamental particle or its decay prod- 
ucts (e.g. Farrar & Biermann, 1998) or quantum-mechanical 
topological defects (such as superconducting strings, e.g. Sigl 
et al. , 1994), there are two major classes of models describ- 
ing the production of UHECRs: 

• Continuous production of UHECRs in the kpc-scale 
shocks associated with the jets of powerful Fanaroff-Riley 
class II radio galaxies (Rachen & Biermann, 1993; Norman, 
Melrose & Achterberg, 1995), or in shocks associated with 
ongoing large-scale structure formation in clusters of galax- 
ies (Norman, Melrose & Achterberg, 1995; Kang, Ryu & 
Jones, 1996; Kang, Rachen & Biermann, 1997); 

• Impulsive production in the same sources respon- 
sible for gamma-ray bursts (Waxman, 1995a; Vietri, 1995; 
Milgrom & Usov, 1995). This possibility has been proposed 
on the basis of a striking coincidence: the typical UHECR 
energy flux above 10 19 eV, and the typical gamma-ray 
flux from gamma-ray bursters in the cosmological scenario 
(Paczynski, 1986) is roughly the same. The recent identifi- 
cation of optical counterparts for several gamma-ray bursts 
(van Paradijs et al. , 1997; Djorgovski et al. , 1997; Metzger 
et al. , 1997) seems to indicate that they indeed occur at 
cosmological distances, and not in an extended halo around 
our own galaxy as had also been proposed. 

A recent review of the various issues involved in the produc- 
tion and propagation of UHECRs can be found in Biermann 
(1995). 

A complication in all these models is the limited vol- 
ume in the local universe, corresponding to a typical ra- 
dius D max ~ 50 Mpc, from which UHRCRs above 10 195 
eV must originate if they are indeed protons or light nuclei, 
and not some exotic particle. Interaction between photons 
of the Cosmic Microwave Background (CMWB) and UHE- 
CRs leads to photo-pion production, e.g. p + 7 — > p + it's. 
This process strongly degrades the energy of particles above 
the Greisen-Zatsepin-Kuz'min cut-off energy E c ~ 10 19 - 5 eV 
(Greisen, 1966; Zatsepin & Kuz'min, 1966) with a typical 
loss length of about 20-30 Mpc above 10 20 eV. 



The observation in 1991 with the Fly's Eye detector of 
a particle with an energy of E ~ 3 x 10 eV (Bird et al. , 

1995) , and more recently of six events above 10 20 eV with 
the Akeno Giant Air Shower Array (Takeda et al. , 1998), all 
well above the GZK cut-off, suggests either relatively close- 
by but unidentified continuous sources (D < 50 Mpc, e.g. 
Elbert & Sommers, 1995), or the scenario where particles 
are produced impulsively and arrive with a strongly energy- 
dependent delay with respect to photons due to scattering 
on weak intergalactic magnetic fields (Waxman, 1995a,b; 
Waxman & Coppi, 1996; Miralda-Escude & Waxman, 1996; 
Waxman & Miralda-Escude, 1996). In the latter case, an 
identification of the source in other wavelength bands is not 
to be expected if the delay is large enough. As we will show, 
for reasonable values of the magnetic field strength in the 
intergalactic medium, the delay of UHECRs with respect to 
photons could be as large as 10 3 — 10 5 years. 

In this paper we will not address the production mecha- 
nism of UHECRs. Rather we will concentrate on the propa- 
gation of these particles in the intergalactic medium (IGM) 
under the influence of weak IGM magnetic fields and their 
interaction with the cosmic microwave background. Scatter- 
ing on IGM magnetic fields leads to diffusion in the flight 
direction of UHECR particles, which induces a time de- 
lay with respect to photons emitted simultaneously at the 
source (e.g. Waxman, 1995; Miralda-Escude & Waxman, 

1996) . The latter effect is especially important in GRB mod- 
els of UHECR production, where a time delay between the 
primary gamma ray photons and 'secondary' UHECRs pro- 
duced concurrently may yield a strong observational con- 
straint on the strength of the IGM magnetic field. Both the 
time delay and rms deflection are strong functions of par- 
ticle energy, t^d <x E~ 2 and a rm s oc E^ 1 , so that losses 
incurred by UHECR particles due to their interaction with 
the CMWB (pion production, and pair production below 
10 19 eV) strongly influence particle propagation. 

In Section 2 we systematically derive the relevant equa- 
tions for the propagation of UHECRs in a weak, random 
intergalactic magnetic field. In Section 3, the precise form 
of the energy losses is considered. We re-examine the treat- 
ment of energy losses of UHECRs and show that a continu- 
ous approximation in terms of the mean energy loss is not a 
good approximation for particles above 10 19 eV. In partic- 
ular we will show that the correct treatment of the Poisson 
statistics of photon-UHECR encounters leads to a signifi- 
cant high-energy tail in the spectrum for sources closer than 
~ 50 Mpc. In Section 4 we present numerical simulations 
of cosmic-ray propagation. The observational consequences 
in terms of the different UHECR production scenarios are 
discussed in Section 5. Conclusions are presented in Section 
6. 



2 UHECR PROPAGATION IN 
INTERGALACTIC FIELDS 

In the tenuous intergalactic medium (IGM) , magnetic fields 
are responsible for the deflection of UHECRs. Direct mea- 
surements of the IGM field strength outside clusters arc 
not available. Observational limits (Kronberg, 1996) can be 
placed on the ordered (large-scale) component (B) IGM and 
the amplitude B T of the random component of the inter- 
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galactic magnetic field. The random component is assumed 
to have a coherence length £ co h — 1 Mpc. These limits are: 



(B) 



IGM 



< 10" 



< 10" 



(3) 



The field strengths within clusters of galaxies can be much 



larger, B ~ 10 



10" 



The typical gyration radius of an ultra-relativistic cos- 
mic ray particle with charge q = Ze and energy E in a 
magnetic field B is 



E 
ZeB 



0.1 



E20 
ZB-t 



Here we have defined 
E 



E20 = 



10 20 eV 



B- 



Gpc . 



B 



10- 



(4) 



(5) 



These estimates imply that the turning angle of UHECR 
particles in the average field (£>} IGM over a distance D will 
be 



D 



r S (E) 



0.53° 



ZD 2 

E20 



(B) 



IGM 



10- 11 G 



(6) 



In this expression we defined D2 = Z)/(100 Mpc). This es- 
timate assumes a uniform field so it is actually an upper 
limit. 

The time delay with respect to ballistic propagation 
follows from expanding the relationship between the path 
length I along the circular orbit of a charged particle with 
gyration radius r g (E) and the linear distance D travelled. 
It equals, for £ <C r g (E), 



^dcla 



D 



24cr 2 (,E) 



1200 



10" 



yr ■ 



(7) 



In what follows, we will assume that the influence of the 
mean field on UHECR propagation can be neglected with 
respect to the action of the random field component. We 
do not believe that a the intergalactic magnetic field will 
be uniform on scales exceeding ~ 10 Mpc, the typical scale 
associated the large-scale structure of our Universe. In that 
case, the contribution of the large-scale field to the delay 
would be measured in terms of a few years around 10 20 eV. 



2.1 Small-angle scattering by random fields 

The equations of motion for an ultra-relativistic particle 
with charge Ze and energy E in a quasi-static magnetic 
field can be written as 



dx 

Til 



dn 

Is 



Ze\B\ 
E 



n x b 



(8) 



Here s = ct is the path length along the orbit, n = 
(ni , ri2 , 713) is the unit vector along the direction of flight 
and 6 = (61 , 62 , &a) is the unit vector along the magnetic 
field. For simplicity we will assume that the magnetic field 
has a uniform amplitude B T and a random (isotropic) direc- 
tion distribution so that 



B r 



i 8 

3 Ul 3 



(9) 



Here the double brackets ((• • ■)) denote an ensemble average. 
This assumption is equivalent to a collection of randomly 
oriented cells with a uniform value of the magnetic field 
strength. The quasi-static treatment of the magnetic field 
corresponds to the requirement 



(n ■ V)6 > 



c dt 



(10) 



which is always satisfied as the bulk (turbulent) velocity in 
the IGM is much less than the speed of light. 

The equation of motion for n can be written in compo- 
nent form as: 



dn; 
~dT 



r s (E) 



(11) 



with r s (E) = EjZeB T and eij k the totally anti-symmetric 
Levi-Cevita tensor. The summation convention is used here 
and below. Magnetic scattering leads to diffusion of the flight 
direction n. The corresponding diffusion coefficient, 



V, 



(An; An.,) 
2As 



(12) 



follows from the Kubo- Taylor formula (Taylor, 1921; see 
also: Sturrock, 1994) as 



Va = 



dse tk in k (0)e ]qr n q (s) (bi(0)b r (x(s))) .(13) 



This diffusion tensor is defined per unit path length, and its 
components have the dimension [length] _1 . If the deflection 
of the particle incurred in one correlation length is small in 
the sense that |An| ~ £ co h/r g (E) 1, the flight direction 
h may be approximated as constant (ballistic motion) and 
rik(0)n q {s) w nfc(0)n 9 (0). This essentially corresponds to 
the well-known quasi-linear approximation of wave-particle 
interactions (e.g. Davidson, 1972). Using Eqn. (^) one has 



ds (6,(0)6r(x( a ))> 



Sir 



(14) 



This integral of the two-point correlation (bi(0)b r (x(s))) for- 
mally defines the coherence length ^ co h of the random field. 
Using the property e ljk e Urn = Oji8 km - 8 jm S k i of the Levi- 
Cevita tensor, the resulting diffusion tensor can be written 
in terms of a scalar diffusion coefficient T>q and a projection 
tensor <5jj — UiUj onto the plane perpendicular to n: 



T>ij = V (8ij — UiUj) 



V = 



3r 2 (£) 



(15) 



It is easily checked that this diffusion tensor preserves the 
unit norm n • n = 1 since n; T>ij rij = 0. 

The expression for the scalar diffusion coefficient T>o can 
be explained in terms of a simple model. A charged particle 
crossing a magnetic cell of size 2£ co j 1 and field strength \B\ ~ 
B r turns through an angle 



8a : 



2ZeB ± £ coh 
E 



(16) 



where B±_ is the component of the field perpendicular to the 
particle velocity. Deflections in different cells add up diffu- 
sively. After encountering N ~ s/2£ co h cells in a distance 
s particles have turned an angle ~ a rm s corresponding to 
Q?ms = (s/24oh) x (((2ZeB±£ coh /E) 2 )). The associated 
angle diffusion coefficient is related to scalar diffusion coef- 
ficient T>o by a 2 ms = 4T>os (see Eqn. B.17), so one finds: 
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T> 



.(«) 



4s 



(17) 



Here we used ((B 2 j_)) = §B r 2 . 

The result ( p^ ) carries over to the more general case ol 
scattering by a broad-band spectrum of random magnetic 
fields. As shown in Appendix A, an isotropic distribution of 
magnetic fields with a power spectrum B(k) as a function of 
wave number, normalised in such a way that 



By 



= <l*l 2 > = 



AkBik) 



(18) 



leads to a scalar diffusion coefficient of the form X>o = 
^coh/3rg(_E) with r g (E) = E/ZeB xxas and an effective corre- 
lation length equal to 



3tt 
8£> rm 



dfc 



B(k) 
k 



(19) 



Here fe m i n ~ 2n/r li (E) is the smallest wavenumber for which 
the approximation of using the unperturbed particle trajec- 
tory still (marginally) applies. 



2.2 Deviation angles and time delays 

In this Section, we consider the deviation angle with respect 
to the line of sight to the source, and the time delay with 
respect to photons which is accumulated by UHECRs due to 
small-angle scattering on magnetic cells in the intergalactic 
medium. Detailed calculations which consider the influence 
of scattering on the flight direction and position of the par- 
ticles, can be found in Appendix B. 

Scattering leads to a deviation angle between the origi- 
nal direction of flight, taken to be along the z-axis, and the 
direction of propagation n. Defining a = cos _1 (n ■ e z ) one 
has after a path length s: 



(a 2 ) (s) ~ AV s = 



(20) 



This differs from the naive estimate, a 2 (s) ~ 2T>os, by a 
factor of two due to the effect of dynamical friction. 

However, the angle between the line of sight from the 
source, f = r/\r\, and the particle flight direction as mea- 
sured by an observer, defined by a' = cos _1 (r ■ n), will be 
less due to the statistical correlation between the particle 
position, r — (x , y , z), and the flight direction n. One finds 
for this angle {a') = and (a') 2 (s) = | Vqs after a path 
length s (Eqn. B22). To zeroth order in T>qs the linear dis- 
tance to the source equals D ~ s (see below). An observer 
at a distance D will therefore measure a spread of UHECR 
arrival directions with respect to the local line of sight from 
the source equal to 



(a') 2 (s^D) = l 



DL 



(21) 



which corresponds to and rms value a rn 
to 

Qrms „ 3.5° (D 2 £ ) 1/2 . 

Here we defined £ = £ coh /(l Mpc) and B_ 9 = B T /(10~ 9 G). 
This estimate neglects the changes in particle energy as it 



(a 1 ) 2 equal 



(22) 



propagates. Energy losses are important, as we will see in 
Section 4. 

The typical angle between the line of sight to the 
source from the observer and the original flight direction, 
6 = cos _1 (r ■ e z ), satisfies (# 2 (s)^ = | T>os, so its rms value 
satisfies 8 lms ~ a rms . 

Due to scattering, the path of an individual UHECR is 
corrugated, leading to a path length increase with respect to 
ballistic (straight line) propagation. The average scattering- 
induced arrival time delay with respect to ballistic propaga- 
tion at t he sp eed of light at distance D, (tdei) = (s — D) / c, 
is (Eqn. |B1E|): 



which is 



D 



(tdel) 



3.1 X 10' 



r s (E) 



V E20 



Di £0 yr . 



(23) 



(24) 



We again assume a constant particle energy. All these rela- 
tions are derived with the implicit assumption that all devi- 
ation angles remain small, i.e. T>os <g; 1. 

The deflection angle a rms and the arrival delay (tdei) 
should be interpreted with some care. Particles originating 
from the same source will only see statistically independent 
realisations of the random magnetic field if they diffuse more 
than one coherence length apart in the direction perpen- 
dicular to the unperturbed orbit. The typical perpendicular 
distance from the original direction of flight at a path length 
s from the source equals 



(25) 



where we use Eqn. ( B17 ). The requirement d± > £ co h deter- 
mines the (energy dependent) decorrelation distance: 



( 3^coh 



1/3 



(|) 2/3 4oh 1/3 r 2 / 3 (S) 



which is 

D dc ^30£^(^- 



2/3 



Mpc 



(26) 



(27) 



Total decorrelation occurs for D 3> -Ddc only. Particles with 
energies such that the source distance is less than Ddc will 
incur a deviation angle of order a rmB , and an arrival de- 
lay with a typical magnitude (tdel)j but the dispersion in 
these quantities associated with a collection of particles from 
the same source is expected to be smaller than for particles 
originating at distances much larger than Ddc, which have 
traversed different magnetic cells. This means that the sim- 
ulation results of Section 4 should be interpreted differently 
for D<Dd c and D S> Dd c - In the first case, the distribu- 
tions are probability distributions for the time delay, devia- 
tion angle etc. averaged over an ensemble of sources. In the 
second case, particles from a single source will be actually 
distributed according to the distributions shown. Given the 
low fluxes above 10 18 ' 5 eV, it is presently not possible to 
measure the resulting dispersion from a single source, and 
only individual arrival directions are measured. 
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2.3 Spatial diffusion 

A description of UHECR propagation in terms of spatial 
diffusion only applies for low energies and/or for large dis- 
tances. 

At low energies, where r g (E) < £ co h, the deflection angle 
in traversing a single cell is already large, 

with E* the characteristic energy for a given field strength 
and field reversal scale: 



(28) 



E, = ZeB T 



9.2 x W 17 ZB- 9 £ eV 



(29) 



Such particles will diffuse through the network of magnetic 
cells, and presumably within cells if the magnetic turbulence 
has a broad band spectrum. To describe their propagation, 
the details of the spectrum of magnetic fluctuations becomes 
important. As a rough estimate, we assume Bohm diffusion 
with a mean-free-path ~ r s (E) = £ co h (E/E*) and spatial 
diffusion coefficient JCb = cr s (E)/3. In that case the finite 
age of the Universe, t ~ H^ 1 , limits the maximum distance 
of sources contributing to the present-day flux at E<E* to 



D diB (E < E, 



rJ 



1/2 / E v 1/2 



55 £,' A h 



1/2,-1/2 



(i) 
(f) 



1/2 



Mpc . 



where h = H /(W0 km s _1 Mpc -1 ). 

This essentially means that the Universe is not trans- 
parent for cosmic rays with E < E* provided it is filled 
uniformly with magnetic cells. This estimate ignores the ef- 
fects of cosmological evolution. These effects are however 
included the simulations presented in Section 4. 

Assuming that UHECRs are indeed of extragalactic ori- 
gin, this implies an exponential cut-off in the UHECR flux 
at Earth below E ~ E*. For particles with E > E t the de- 
viation angle accumulated by traversing many cells becomes 
large (see Eqn. B22/23) at a distance of order 



X(E) 



1 _ H(E) 
2T>o 2£ co h 



which is 



HE) = 15 (jj^j *o 1 Gpc 



(30) 



(31) 



As shown in Appendix B, the propagation of UHECRs on 
these scales can be treated as spati al diffusion with a diffu- 
sion coefficient equal to (Eqn. B25 ) 



K : 



D 2 

-^=cX(E) 



(32) 



The finite age of the universe, t ~ H , gives a maximum 
distance of UHECR sources contributing to the present-day 
flux: 



D dia (E) 



(33) 



The rapid energy losses due to photo-pion production above 
10 19 ' 5 /(1 + z) eV, followed by a much slower loss at lower 
energies, means that effectively one can put E < 10 19 ' 5 eV 



in these estimates. So for E > E„, sources contributing to 
the UHECR flux must be closer than 



Aiiff(10 eV) ~ 600 £\ 



-1/2,-1/2 



(ZB-g)- 1 Mpc . 



(34) 



2.4 Effects of large-scale structure 

So far, we have treated the effects of the intergalactic mag- 
netic field in a simple model where the field consists of ran- 
domly oriented cells which are distributed more-or-less ho- 
mogeneously. Although true observational evidence is lack- 
ing, one might argue that magnetic fields follow the dis- 
tribution of luminous matter and are concentrated mainly 
in the filamentary large-scale structure of the universe (e.g. 
de Lapparent et al. , 1986). A concentration of strong mag- 
netic fields in the filaments observed in the large-scale struc- 
ture could be expected, in particular if fields are amplified 
considerably in the turbulence associated with the forma- 
tion of galaxies, as proposed by Kulsrud et al. (1997). Such 
structure in the field distribution is important for the de- 
lay incurred by secondary photons produced by the decay 
of photo-produced neutral pions (Waxman & Coppi, 1996): 
in the impulsive scenario, the delay of secondary photons 
will mirror the delay incurred by the UHECRs at their pro- 
duction site with respect to the primary photons associated 
with the explosive event. 

If the IGM magnetic field distribution is strongly struc- 
tured on scales of > 10 Mpc, the observational upper limit 
for the typical field strength, B T < 10~ 9 G, does not apply. 
Random fields on Mpc scales, concentrated in cluster or su- 
perclusters, could be as strong as 10 -7 — 10 -6 G. In that case, 
a totally different model of UHECR propagation should 
be considered, where particles at lower energy are trapped 
within clusters in a manner reminiscent of the 'Leaky Box' 
models often used in the description of the propagation and 
confinement of Galactic cosmic rays (Rachen, 1998 private 
communication). Such a description would extend to parti- 
cles with energy E < 10 20 Z (B r /0.1 ^G) (i/Mpc) eV, where 
L is the size of the cluster. 



3 UHECR ENERGY LOSSES 

UHECRs with energies above 10 18 5 eV lose energy pri- 
marily through reactions with photons of the cosmic mi- 
crowave background (e.g. Berezinskii et al. , 1990 and refer- 
ences therein). Because of the strong energy dependence of 
scattering on IGM magnetic fields, T>o oc E~ 2 , this energy 
evolution has to be considered in detail when describing the 
intergalactic propagation of UHECRs. We parameterise the 
mean energy loss per unit path length by the energy loss 
length 1(E) defined by 



£{E) = 



E ~ds 



(35) 



At the lowest energies under consideration in the simulations 
we will describe in the next Section, 10 17 < E < 10 18 eV, the 
dominant loss mechanism is (adiabatic) expansion losses in 
the Hubble flow. The loss length equals the Hubble length, 



1 Gpc . 



(36) 
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For protons in the energy range 10 18 ' 5 < E < 10 19 ' 5 eV, the 
dominant loss mechanism is photon pair production: p+7 — > 
p+e + + e~ (Blumenthal, 1970). The energy loss length £ P (E) 
for this process reaches its minimum value in the energy 
range 10 185 < E < 10 19,5 eV (see the figure below and also 
fig. 3 of Yoshida & Teshima, 1993): 

e p (E) ~ 2 Gpc . (37) 

At higher energies, losses are dominated by photo-pion pro- 
duction: p + 7 — > p + 7r's. As more and more CMWB pho- 
tons in the Wien tail of the blackbody spectrum satisfy 
the threshold condition for pion production at increasing 
proton energies, the associated loss length £ W (E) decreases 
exponentially with increasing energy. In the energy range 
10 19,5 < E < 10 20 - 5 eV it scales as (see the discussion be- 
low) 

M£) ^ 4.8 (^-] 2 e E ^ /E Mpc , (38) 

where E th ~ (m 7r m p c i ) / (2k h TcMWB) — 3 x 10 20 eV is the 
threshold energy. For E > 10 20 ' 5 eV pion production on the 
CMWB saturates at a constant loss length, as essentially 
all CMWB photons exceed the threshold in the proton rest 
frame: 

e„(E) ~ 10 Mpc. (39) 

The energy loss length in the local universe (z = 0) resulting 
from these combined processes is shown in Figure 1. 

In order to describe the UHECR spectrum received at 
Earth, one has to calculate the modification of the spectrum 
at the source due to pion-production on the CMWB. A num- 
ber of pioneering studies have been done in the past (e.g. 
Berezinskii, Grigor'eva & Zatsepin, 1975; Hill & Schramm, 
1985; Berezinskii & Grigor'eva, 1988; Rachen & Biermann, 
1993; Yoshida & Teshima, 1994 and Aharonian & Cronin, 
1994). 

Much of the discussion has focussed on the approxima- 
tion to be used to describe the spectrum near the Greisen- 
Zatsepin-Kuz'min (GZK) cut-off in the spectrum (Greisen, 
1966; Zatsepin & Kuz'min, 1966), which results from the 
exponential dependence of the photo-pion production loss 
length (Eqn. K8() . This cut-off occurs at an energy E c ~ 
10 19 ' 5 eV. Particles injected above the cut-off are expected 
to 'pile up' close to E c , forming a bump in the spectrum. 
The height and width of this bump is determined by the 
details of the pion production process. 

Early work (e.g. Berezinskii & Grigor'eva, 1988) as- 
sumed that the effect of pion losses could be described by a 
simple continuous energy loss approximation, i.e. di?/ds = 
-E/l{E). It was pointed out by Hill & Schramm (1985) that 
there is an intrinsic spread in energy which is not described 
by this simple approximation. This spread results from two 
effects: 

(i) The kinematics of the photo-pion production, which 
leads to a spread in the energy of particles in the observer's 
frame due to the spread of pion production angles in the 
center-of- momentum frame (CMF). 

(ii) The Poisson noise in the number of photons encoun- 
tered for a given path length s. 

These effects modify the result of the continuous loss ap- 
proximation, spreading the bump over a larger energy range 



energy loss length 




- H„ = 75 km s 1 Mpc 
! q = 0.5, z = 



10 17 10 18 10 19 10 30 io 21 10 22 

E(eV) 

Figure 1. The loss length 1(E) (thick solid curve) as a function 
of energy. It is the result of losses due to pion production (in, 
dotted curve), pair losses (£ p , dashed curve) and expansion losses 
in the Hubble flow (In = c/Ho, horizontal solid line) in a flat 
Universe (go = 1). 



and making it less prominent. In the next section, we will 
consider a simple model for the evolution of the UHECR 
spectrum in transit from source to observer due to these 
effects. 



3.1 Kinematics of photo-pion production 

Because of its importance in determining the GZK cut-off 
energy, and the 'survival probability' of UHECRs above 
the cut-off, we consider photo-pion production in some- 
what more detail. A more complete account can be found 
in Berezinskii et al. (1975) or in Mannheim & Biermann 
(1989). All expressions are for protons (Z — A = 1), but 
can be readily generalised to other nuclei. We use units with 
c = 1 in this Section. 

The interaction between a proton with energy E p and a 
photon with energy e in the observer's frame is most conve- 
niently described in terms of the invariant total energy E t 
in the center of momentum frame (CMF), which moves with 
Lorentz factor 7cmf = (E p + e)/E t ~ E p /E t 2> 1, and the 
photon energy eo in the proton rest frame (PRF). They are 
related by 

E t = \J m p 2 + 2m p e , eo = 7p £ (1 - Pp cos 8) . (40) 
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Here j p = — /3 P 2 is the Lorentz factor of the proton, 

and 6 the angle between the photon propagation direction 
and the proton momentum in the lab frame. A head-on col- 
lision between proton and photon corresponds to 8 = n. 

Energy and momentum conservation in the CMF, fol- 
lowed by a Lorentz transformation back to the observers 
frame, yield the following expression for the final proton en- 
ergy: 

E p l ~E P [l-Kp + iC cos crjf ] . (41) 

Here an is the angle between the momentum vector of the 
incoming and final proton in the CMF. We have put /3cmf ~ 
1, which is an excellent approximation at these energies, and 
defined a mean elasticity 



2E t 



(42) 



The quantity K is a measure of the spread of final energies 
around the mean, 



K 



(E* - m+ 2 )(E t 2 - m^) 



2E t 



(43) 



with m± = m p ± m^. If one assumes that the CMF produc- 
tion angle an is distributed isotropically, the mean energy 
loss per interaction in the observer's frame equals AE P — 
—K P E P and the energy dispersion is A_B r ms — (K/yo) E p . 

The spread in the final proton energy in the observer's 
frame increases strongly with CMF energy. Near reaction 
threshold, which corresponds to zero momentum for the final 
proton and the production of a single pion at rest so that 
Et — m p + m.7r, one has K p ~ 0.13 and K « 0, and the 
spread in the final proton energies in the observer's frame 
is small. Much above threshold, where Et 3> m+ and K p « 
K « ~, one has E p w | E p (1 + coscif), so the spread in 
final energies is large. 

The interaction rate for photo-pion production in the 
observer's frame, 1Z P1 , can be expressed in terms of the reac- 
tion cross-section <r P7 (eo) in the PRF by using the Lorentz- 
invariance of the photon occupation number in phase space 
and the relativistic aberration formulas which, in the limit 
7 P ^> 1, lead to strong relativistic beaming of the CMWB 
photons in the PRF, so that all photons have cos So — — 1. 
The relevant expressions can be found in Blumenthal (1970). 

This leads to an expression for the reaction rate derived 
by Berezinskii et al. (1975), but with corrected integration 
limits corresponding to the kinematically allowed values: 



de 



n(e)/j_ 
t 2 ^7p 2 



de ^ 7 (eo)eo ) .(44) 

th/27p 

Here e t h = rn n (l + m n /2m p ) ~ 145 MeV is the photon 
threshold energy in the PRF. Berezinskii et al. (1975) incor- 
rectly replace the upper limit on the integration over eo for 
fixed e by infinity. For a thermal (Planckian) photon distri- 
bution in the observer's frame with temperature T one has, 
reinstating c: 



n(e) 

c 2 



1 



e/k b T 



The PRF pion production cross-section <7 P7 (eo) 



(45) 



near 
o 



reaction near s 



E t 



1.6 GeV 2 



with 0" P7 (s m ) 



0.5 mb. Therefore, not too far from threshold, we can ap- 
proximate the cross-section by 



y (s) ~ As a P7 (s m ) <5(s - s m ) 



(46) 



in the integral over eo or s. Here As ~ m p e t h — 0.15 GeV 
is the width of the resonance. The value of s m corresponds 
to e ~ 0.3 GeV » 2e th . 

Using the above approximation for the cross-section, the 
reaction rate for E p <C E t h typically equals 



Eth 

~eZ 



~ E th/ e p 



(47) 



with a m — a PJ (sm), and E t h = (m 7r m p c 4 )/(2fcbT). Putting 
Eth 2 = s m in Eqns. @ and @ yields K p ~ 0.2 and 
K ~ 0.165 in this limit. 

At energies E p 3> -Eth> Et S> m p + m w , the cross-section 
summed over all reactions becomes approximately constant, 
a 7T ~'(s) ~ 0.16 mb = ao- The reaction rate 7?. P7 (-B P ) then 



approaches the constant rate 



R P -y(E p > Eth) ^ 0.244 (^-) i .vt , . 



(48) 



while the mean fractional energy change and energy disper- 
sion satisfy K p ~ K ~ 0.5. The energy loss length 1(E) 
follows from these quantities as 



£(E) = c/ (1l pl (E) K p ) 



3.2 Effects of Poisson noise 



(49) 



The mean number of pion-producing photons encountered 
by UHECR particles in a path length As = cAt equals 



(n P h) (As) 



As 



K P £(E) 



(50) 



threshold is dominated by a resonance in the p+7 — > p+iv 



The actual number of photons encountered by an individual 
particle is subject to Poisson statistics. This means that the 
variance in the number of encounters equals \n = (nph) 1 ^ 2 , 
while the probability of a particle not encountering any pho- 
tons decays as P(n p h = 0) = e~( Uph ). The effect of Poisson 
statistics is especially important for protons at high energies, 
E p > 10 18 5 eV, where K v l(E) ~ 20 Mpc. For protons origi- 
nating from sources closer than D ~ 20 — 50 Mpc, the effect 
of Poisson statistics on the energy loss should be clearly vis- 
ible in the arrival spectrum. There is a significant fraction 
of particles which have interacted with no (or only a few) 
photons, leading to a high-energy tail in the spectrum which 
more or less reflects the original source spectrum, because 
the loss length 1(E) is roughly independent of energy, but 
at an amplitude reduced by a factor ~ e ~{ n ph)( D \ 



4 SIMULATIONS 

We have constructed a simulation code for the propagation 
of protons in intergalactic space, taking into account the 
scattering on IGM magnetic fields and the energy losses suf- 
fered as a result of interactions with the CMWB. 

The propagation of the UHECRs over a path length As 
is calculated with the simple explicit scheme 
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x(s + As) = x(s) + n(s) As 



(51) 



inelasticity 



The unit vector n along the direction ol propagation under- 
goes diffusion, which is described by the the explicit numer- 
ical integration of a stochastic differential equation (SDE). 
The stochastic change in the direction of flight is modelled 
by 



An at = V2V As +£ 2 e 2 ) 



(52) 



Here £i and £2 are two independent unit Wiener processes, 
drawn at each step from a Gaussian distribution with unit 
dispersion so that ($1,2) = 0, (£1,2) = 1- The two (arbitrary 
but mutually orthogonal) unit vectors ei ]2 lie in a plane 
perpendicular to n so that An st • n = 0. Note that n = 
(ni , n 2 , ns) and ei :2 are the eigenvectors (with eigenvalues 
and 1 respectively) of the diffusion tensor T>ij — T>o Pij, 
with Pij = (Sij — mrij), which motivates this approach. This 
prescription is equivalent to the stochastic term in Eqn. ( Bl ) 
of Appendix B, with the correspondence 



P(n) ■ AW — > V2V As (£iei + £ 2 e 2 ) . 
The new propagation direction then is calculated as 



n(s + As) = \ 1 - I An 



+ An st 



(53) 



(54) 



This preserves the norm n • n = 1 and is an excellent ap- 
proximation for the diffusion process provided |An st | <C 1. 
When many independent realizations of these SDEs are sim- 
ulated, the resulting distribution of unit vectors corresponds 
to a solution of the diffusion equation for n. Details of similar 
astrophysical applications of SDEs can be found in Achter- 
berg & Kriills (1992). 

The change in UHECR energy is calculated as the com- 
bined effect of expansion losses in the the Hubble flow, 
pair production and photo-pion production losses. Pair pro- 
duction is treated using the continuous loss approximation, 
(dE/ds) p = —E/£ p , which is an excellent approximation for 
this process. We use the interpolation formula of Stepney & 
Guilbert (1983) for the cross-section as a function of proton 
energy, and the numerical fit of Chodorowski, Zdziarski & 
Sikora (1992) for the inelasticity of the reaction. The inte- 
gration over the CMWB photon distribution in the proton 
rest frame proceeds along the lines elaborated in Section 
3. The pair production loss length £ P {E) so obtained (see 
Figure 1) closely follows the result of Blumenthal (1970). 

Pion production must be treated in more detail. The 
number of pion-producing photons encountered in a path- 
length increment As is drawn from a Poisson distribu- 
tion with mean (n p h) = As/£ pl . The interaction length 
l pl = c/lZp-t i s calculated according to (with threshold en- 
ergy B th ~ 3 x 10 20 eV): 



li.gf'-J^ e Bth/B (E<0.2E t 

- Eth J 



Vy(Mpc) 



(55) 



4.8 



(E > 0.2£ t h) 



Here we have written E rather than E p for the proton en- 
ergy. The energy loss per interaction is calculated using a 
simple interpolation formula for the mean inelasticity: 



P V E th + E ) 



(56) 
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Figure 2. The mean inelasticity (K) = K p (solid line) and the 
upper and lower limit (dashed lines), corresponding to (K) ± K, 
due to the variance in production angle CT;f in the PRF. At 
given energy, the inelasticity is uniformly distributed between the 
dashed lines. 



The pion production angle in the center of momentum frame 
is drawn from a uniform distribution — 1 < coscrif < 1. The 
energy change is calculated according to Eqn. (^), with K 



calculated by re-expressing E t 
exact kinematic relations for K B 



in terms of K p using the 
and K: 



K = ^(K P + K+)(K P -K-) , (57) 

with K± = m^/(m p =pm^). This allows the use of a single in- 
terpolation for K P (E) while preserving the correct behaviour 
of K at both low and high energies. This is illustrated in 
Figure 2. 

Finally, at large source distances cosmological effects 
are taken into account. Due to the higher CMWB photon 
density and temperature at a distance corresponding to a 
redshift z, n p h oc (1 + z) 3 and Tcmwb oc (1 + z), the energy 
loss length due to interaction with CMBW photons at red- 
shift z can be obtained from the value in the local universe 
[z = 0) by using the scaling law 



\{E, 



e{(l + z)E,z = 0) 
(1 + z) 3 



(58) 



This is a simple consequence of the fact that the loss length 
is inversely proportional to the photon number density and 
depends on energy E and the lab frame threshold energy 
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energy distribution 



D = 10, 25, 50, 100 Mpc 
E in] . = 10 al eV 

= 10" 10 G, <„„,, = 1 Mpc 



« 'o 

X " H 

n 




Figure 3. Arrival spectra at Earth for mono-energetic injection 
at an energy of 10 21 eV at the source, for four source distances D. 
The spectra have been multiplied by an arbitrary factor oc D 2 . 
For increasing distance, the maximum of the spectrum shifts to 
lower energy. For 10 Mpc there is still a significant fraction of 
particles at the injection energy, giving rise to a distinct spike. 
This spike has almost completely vanished for D = 50 Mpc. The 
spectrum above 10 19 ' 5 eV, where the loss length t(E) decreases 
by almost two orders of magnitude, decreases rapidly for larger 
distances. 



E t h oc Tcmwb oc (1 + z) 1 of the interaction through the 
ratio E/Eth oc 1 + z. 

The random magnetic field strength in the general in- 
tergalactic medium scales as B T oc (1 + z) 2 as a result of 
flux freezing in an expanding universe (e.g. Subramanian & 
Barrow, 1998). The universal expansion stretches the cor- 
relation length so that £ co h(z) = £ co h(z = 0)/(l + z). As 
a result of these effects the magnitude of the flight direc- 
tion diffusion rate scales as T>o oc B 2 £ co ^ oc (1 + z) 3 . In our 
simulations we assume a flat, matter-dominated universe in 
which the Hubble length scales as in = (c/Hq) x (1 + z)~ 3//2 . 



4.1 Simulation results 

To illustrate the effects of UHECR scattering most clearly, 
we consider the case of mono-energetic injection first. Parti- 
cles are injected with an energy Ei n j = 10 21 eV. Arrival spec- 
tra, delay times and angular distributions are calculated for 
different source distances. Figure 3 shows the arrival spectra. 
From this calculation it is obvious that the effect of 



Poisson statistics remains strong up to distances of about 50 
Mpc: a tail of high-energy particles up to the injection en- 
ergy remains. At larger distances, the probability of not en- 
countering any CMBW photons becomes vanishingly small. 
Energy losses proceed much more slowly once particles reach 
an energy E < 10 19 ' 5 eV, as the loss length becomes larger 
than 1 Gpc. Therefore particles tend to accumulate around 
10 19 ' 5 eV, drifting slowly to lower energies. 

Figure 4 shows the distribution of the delay time with 
respect to simultaneously produced photons, for a number 
of different source distances. At small distances they are 
centered around a value close to the predicted value of idei 
for the injection energy of 10 21 eV. For those distances larger 
than 25 Mpc, the delay distribution has its maximum close 
to the value of idei associated with the typical arrival energy, 
£ arr ~ 10 19 - 5 eV. 

Figure 5 shows the cumulative angular distribution of 
the UHECRs for the arrival angle a' with respect to the 
line of sight to the source. The deflection process consid- 
ered here will lead for constant particle energy to an angu- 
lar distribution of the form dN/dQ oc exp(— a' 2 /a 2 ms ), with 
df2 = 27rsina' da'. For small angles, a' and a rms <C 1, this 
distribution scales as 



dAf 
din a' 



(59) 



which has a maximum at a' — a rms . The simulated distribu- 
tions shown in Figure 5, calculated with energy losses taken 
into account, closely resemble model distribution (^). 

Next we consider the case of power-law injection, where 
particles are produced at the source according to 



dN inj {E) = n ini (E) dE = k E~ s dE . 



(60) 



Figure 6 shows the arrival spectrum at Earth as expected for 
a steady source. The injection spectrum was n(E) oc E~ 2 
between 10 17 and 10 21 eV. The delay incurred in transit 
has a negligible influence on the shape of the spectrum as 
long as ctdci <§C D. Therefore the arrival spectra calculated 
for B rms = 10~ n G and £ coh = 1 Mpc are representa- 
tive for all simulations we performed at energies well above 



10 ZB-g£o eV. The spectra have been multiplied 



with a factor oc D 2 E 2 so that the injection spectrum in all 
cases corresponds to a horizontal line. For energies below 
10 19 eV, the effect of losses is negligible at these source dis- 
tances. Particles do accumulate at the GZK cut-off energy, 
but the height of the 'bump' is not very dramatic due to 
the natural spread in arrival energies, AE/E ~ 0(1) for a 
given injection energy due to the variance in the photo-pion 
production kinematics. The high-energy tail above the GZK 
cut-off is very pronounced at D = 10 and at D — 25 Mpc, 
still quite visible in the 50 Mpc curve, but almost absent at a 
source distance D = 100 Mpc. Our results are qualitatively 
similar to those obtained by Yoshida & Teshima (1993) in 
the corresponding distance range. 

Figure 7 shows the delay time idei and the dispersion 
Atdei (shown with error bars as idei ± Aidei) for particles 
originating at D — 25 Mpc as a function of the arrival en- 
ergy _E arr at Earth. The delay time has been scaled by a 
factor oc (E aTT /10 20 eV) 2 to scale away the dominant energy 
dependence predicted by Eqn. (p3|). Particles well below the 
GZK cut-off (E < 10 19 eV) have lost little energy, and there- 
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arrival delay distribution deflection angle distribution 




Figure 4. Delay time distribution at Earth for particles mono- 
energetically injected at E = fO 21 eV for four source distances 
ranging from 10 to 100 Mpc. The distribution for the smallest 
distance is the left-most curve, for the largest distance the right- 
most curve. The maximum in the 10 Mpc curve coincides with 
the theoretical result of Eqn. ^Jfor that distance at the injection 
energy, ^^(lO 21 eV) ~ 0.3 year. The maximum for the 100 Mpc 
curve coincides roughly with the theoretical result for that dis- 
tance and the energy where most particles reside (see figure 3): 
E~2x 10 19 eV, yielding t del ~ 10 5 year 



fore fall almost on the horizontal line predicted by Eqn. (24) 
for Z = 1 (protons): 

t dcl (E, D) V {E 20 /B- 9 D 2 ) 2 ~ 0.31 Myr . (61) 

At the lowest energies considered here, E ~ 10 17 eV, there 
is an increase in the delay time because the small deflection 
angle approximation starts to break down: a rms ~ 20° at 
these energies. Well above the GZK cut-off the delay time is 
roughly half the expected value based on the arrival energy 
due to the fact that the mean energy in transit is larger that 
the arrival energy. For energies near the injection energy, 
E ~ 10 2 eV, the delay time once again approaches the 
theoretical value as these are particles that happen not to 
have interacted significantly with CMWB photons, and have 
therefore propagated with an almost constant energy. 

Figure 8 shows the dispersion in delay times, Cdelj as a 
function of the delay time for the same arrival energy bins 
as in Figure 7, with the dominant energy dependence scaled 
away. From this figure one can conclude that, at least in the 



Figure 5. Angular distribution at Earth of the arrival angle 
Q arr = a! with respect to the line-of-sight to the source, for 
particles monochromatically injected at E = 10 21 eV for dis- 
tances ranging from 10 to 100 Mpc. The maximum of these 
curves is roughly at a rms (Eqn. B2 ). For a distance of 10 Mpc 
and E ~ _B; n j = 10 21 eV, theory predicts ct rm s — 0.012 degrees. 
For a distance of 100 Mpc and an energy of E ~ 2 X 10 19 eV, one 
expects arms — 2 degrees. 



regime considered here, the often-used approximation (e.g. 
Miralda-Escude & Waxman, 1996) 

a del (E, D) = yj((t dc i- < t dci >) 2 ) ~ t dc i(E, D) (62) 

is a fair one. We find roughly a dc \ w 0.6 (idei), but this rela- 
tion is approximate, and breaks down when the delays (and 
deflection angles) become large (not shown). 

Figure 9 shows the angle a' = cos (n • r) with respect 
to the line of sight to the source as a function of arrival 
energy. The dominant behaviour, a' oc -E -1 , has been scaled 
away. If particles did not lose energy, the arrival angle would 
satisfy (Eqn. 22) 

Q ' X (fe) ( D ^r 1/2 -3.5° . (63) 

The actual value (using E = -E arr ) is somewhat less as a 
result of the energy losses. 
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flux/energy distribution 



mean delay 




E (eV) 

Figure 6. Energy distribution n(E) X E 2 , received at Earth for a 
continuously radiating source. Particles are injected at the source 
with n inj (_E) <x E~ 2 between 10 17 < E < 10 21 eV. Spectra are 
calculated for a source distance D equal to 10, 25, 50 and 100 Mpc, 
and multiplied by an arbitrary factor oc D 2 . Only the portion of 
the spectrum between 10 18 ' 5 and 10 21 ' 5 eV is shown. This figure 
should be compared with figure 3 for monochromatic injection. 
There is a clear 'bump' near the GZK cut-off at E ~ 10 19 ' 5 eV. 
Below the GZK cut-off, particles lose energy very slowly (£(E) > 1 
Gpc), and the spectrum is almost independent of source distance. 




E (eV) 

Figure 7. Delay time as a function of arrival energy. The dis- 
persion in the delay time is shown as error bars. The delay has 
been scaled such that propagation at constant energy E would re- 
sult in a horizontal line at tde\ x -^e c v / (-^MpcBnG) 2 — 0.31 Myr, 
cf. Eqn. 24. One sees that for the particles near -E m ax = 10 215 
eV, which happen to have lost little energy, and for the particles 
in the range E < 10 19 ' 5 eV which, at this source distance, have 
lost very little energy, this theoretical value is approximated. For 
10 20 eV < E < 10 21 eV the influence of energy losses reduces 
*del with respect to the value predicted on the basis of the arrival 
energy E arr . 



5 BURSTING SOURCES AND THE GZK 
CUT-OFF 

Observations in the northern hemisphere show the arrival 
direction of events above 4 x 10 19 eV to be distributed over 
the whole sky, but with some evidence for clustering of ar- 
rival directions around the supergalactic plane (Stanev et 
al., 1995; Hayashida et al., 1996). A similar analysis of the 
less extensive data set from the southern hemisphere (Kew- 
ley et al, 1996) does not show evidence for this clustering. 

Due to photo-pion production on the CMWB (Greisen, 
1966; Zatsepin & Kuz'min, 1966) it is likely that particles 
with energy above 10 eV originate from a volume with 
radius D < D m ax — 50 Mpc. Deflection by random inter- 
galactic magnetic fields over a distance D < D max is not 
capable of deflecting UHECRs in this energy range by more 
than 5—10 degrees, given the observational limit on the field 
amplitude B IU1S . 

This means that the observed UHECRs in the energy 



range above 10 19 ' 5 eV (some ~ 100 events, with about 7 
events above 10 20 eV) must originate from multiple sources. 
Solely on the basis of simple energy arguments (Rachen 
& Biermann, 1993; Norman, Melrose & Achterberg, 1995), 
only the most luminous radio galaxies (FRII sources) or 
shocks associated with ongoing large-scale (~ Mpc) struc- 
ture formation can continuously produce particles in this 
energy range. Production sites in the case of radio galax- 
ies are the so-called hot spots, where the impact of large- 
scale collimated jets on the intergalactic gas creates strong, 
large (~ kpc scale) shocks. The typical space density of 
Fanaroff-Riley II radio galaxies in the local universe is 
small, n s ~ 10~ 7 Mpc -3 . The corresponding distance scale, 
D s = (47rn s /3) _1//3 ~ 100 Mpc, is large compared with the 
typical loss length for particles above 10 19 ' 5 eV, £{E) ~ 10 
Mpc. If there is still ongoing structure formation in the lo- 
cal universe, the sky filling factor of the filamentary shocks 
associated with this process would be small. 

It is therefore not surprising that no clear candidate 
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Figure 8. Delay time dispersion <T de i as a function of mean delay 
time (tdel) f° r a source at a distance of 25 Mpc in a 0.01 nG field. 
The dominant energy dependence <x E~ 2 has been scaled away. 
The normalization of the delay time and delay dispersion is as in 
figure 7. At this distance and field strength, where the small-angle 
approximation for the deflection process still applies, one can fit 
the delay dispersion roughly with a linear relation of the form 
""del ~ 0.6 <t dc i>. 



sources for UHECR events have been identified (Elbert & 
Sommers, 1995). Even if one broadens the candidate sources 
to include all active galaxies, n s ~ 10" 5 Mpc" 3 , D s ~ 20 
Mpc, one is hard-pressed to find a sufficient number of can- 
didate sources within the local volume bounded by D ma x, 
the total number of sources scaling as 



N S (D < D„ 



(64) 



A consequence of this is that the UHECR spectrum 
at Earth is always dominated by a few close-by sources in 
the case of continuous production. This is illustrated in fig- 
ures 10 and 11. These show two representative cases of the 
UHECR energy distribution, obtained by a discrete distri- 
bution of sources within a distance of D ln — 250 Mpc, each 
source distance drawn by Monte Carlo simulation from a 
distribution satisfying Eqn. ([ii]). The figure also shows a 
calculation of the background contribution due to sources 
at distances from 250 Mpc to 2.5 Gpc. This background has 
been calculated using a continuous approximation for the 
source density. The source density assumed equals n B = 10~ 6 
Mpc -3 , so that D B ~ 70 Mpc. Particles originating at dis- 



r (eV) 



Figure 9. The arrival angle with respect to the line of sight to 
the source as a function of the arrival energy E aTT . The dominant 
parameter dependence oc (B^J Dl co ^)/ E has been scaled away 
so that particles propagating with constant energy would lie on 
a horizontal line a 1 X (-E20 / -B-gV D2<?o) = 3.5°. Particles have 
been injected with energies uniformly distributed over the range 
10 18 < E < 10 21 eV at a distance of 25 Mpc, and propagate 
through a random field with B rm s = 10 -11 G and £ co h = 1 Mpc. 
The dispersion in the arrival angle is shown as bars. Particles with 
E < 10 19 eV which lose very little energy for this source distance 
follow the theoretical result for constant energy propagation. 



tances larger than 2.5 Gpc will not reach Earth within a 
Hubble time due to the large accumulated delay. These fig- 
ures show that the high-energy part of the spectrum above 
the GZK cut-off at 10 19 ' 5 eV is almost completely dominated 
by the single closest source. 

This simple argument lends some additional credence 
to the gamma-ray burster hypothesis of Waxman (1995a), 
Vietri (1995) and Milgrom & Usov (1995), where UHECRs 
are produced in bursts. 

An important difference with the case where UHE- 
CRs are produced continuously is the fact that in this case 
the number of sources contributing to the UHECR flux at 
Earth at any given time depends critically on the time delay 
tde\(E, D) incurred during propagation: each source remains 
visible for a time At ~ idci as a result of the spread in ar- 
rival times. At any one time, one sees those sources that 
have burst within a time interval At ~ t dc i, so that 



dA?s ~ 47rD 2 Qgrb t de i(E, D) dD 



(65) 
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flux from continuous sources 




Figure 10. The UHECR spectrum due to continuous sources 
with a source density of 10~ 6 Mpc - . It is assumed that all 
sources are identical in injected flux, and that they are distributed 
homogeneously in space so that N(< D) oc D 3 . A Monte-Carlo 
realization with discrete sources is used for D < D m = 250 Mpc. 
The background of distant sources is calculated using a continu- 
ous distribution of source distances between 250 Mpc and 2.5 Gpc 
for the simulated UHECRs. The injection spectrum is a power- 
law, N(E) = k E~ 2 - 5 . The dominant power-law behaviour has 
been scaled away. The thickest line is the total spectrum. The 
two thinner lines are the contribution of the sources within 250 
Mpc and of the background due to more distant sources, which 
only contrinutes below the GZK cut-off at 10 19 ' 5 eV. The three 
thinnest curves are the spectra of the three closest sources in this 
realization of the source distribution. It is seen that the closest 
source completely dominates the spectrum above the GZK cut-off. 
At energies below 10 17 ' 5 eV the closest source again dominates 
the flux as particles from far-away sources do not reach the ob- 
server within a Hubble time due to the slow spatial diffusion of 
particles in this energy range. 



sources contribute to the flux at energy E from a spherical 
shell of thickness dD at distance D. Here Qgrb is the typical 
gamma ray burst rate per unit volume in the cosmological 
scenario. The typical value of Qgrb in a non-evolving sce- 
nario, where GRBs are distributed homogeneously within a 
volume corresponding to a redshift z < 1, is of order 

Qgrb ^ 1(T 8 ( — t) Mpc" 3 yr- 1 . (66) 

VGRB ^kms-iMpc- 1 / y y ' 

The total energy yield per burst is typically £grb — 



flux from continuous sources 




Figure 11. As in figure 10: the UHECR spectrum due to con- 
tinuous sources with the same source density, n s = 10 — 6 Mpc -3 , 
but with a different Monte-Carlo realization of the source distri- 
bution within a distance D; n = 250 Mpc. In this realization there 
are a few more close sources. This distribution is compatible with 
existing observations, whereas the spectrum shown in figure 10 is 
not. 

10 51 - 10 52 erg, assuming a conversion efficiency of 10 % 
into gamma rays. 

However, Wijers et al. (1998) have pointed out that 
most GRB formation scenarios, such as neutron star mergers 
(e.g. Paczyhski, 1986; Narayan, Paczyhski & Piran, 1992) 
or hypernovae (Paczyhski, 1997) involve the end-product of 
massive star evolution. One should therefore expect that 
the GRB rate closely follows, on cosmological timescales, 
the star formation rate. The latter seems to peak at a larger 
redshift (z > 2) . This evolving GRB scenario implies a larger 
distance scale (z > 2) and total energy yield (Egrb ~ 10 53 
erg) for a typical burst, leading to a correspondingly lower 
burst rate per unit volume in the local (z ~ 0) universe. On 
the basis of the observed rate of GRBs, Wijers et al. estimate 
a local burst rate equal to: 

( ^^^ j^Pc-yr-'. (67) 

The local burst rate is the relevant quantity for the produc- 
tion of UHECRs, which must originate within a distance 

D max • 

The main uncertainty in these estimates is the unknown 
beaming factor of the gamma radiation, which could increase 
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the burst rate significantly above the one derived (from ob- 
servations) assuming isotropic emission. In the hypernova 
scenario in particular, where the energy is generated in a 
short-lived torus around a stellar core that has collapsed into 
a black hole, beaming could be significant. Finally, given the 
limited number of events with X-ray, optical and/or radio 
counterparts for which a distance can be derived, one can 
at this time not exclude the possibility that both the hyper- 
nova and binary merger mechanism produce a population of 
GRBs. 

One can define an (energy-dependent) critical distance 
D b (E) which bounds the typical volume within which one 
expects to find a single source, N B (D < D b ) — 1 for given 
energy fi: 



j-D b 

a ' 



4ttL> 2 Qqrb t dcl (E, D) AD = 1 



(68) 



Using idol oc fi 2 it is easily checked that this corresponds to 
(Miralda-Escude & Waxman, 1996) 



— <2grb D b tdei{E, D b ) = 1 
Using Eqn. (24) one finds: 

2/5 



D h (E) ~ 17 



E20 
ZB-< 



-l/B n -l/5 



Q-' J Mpc , 



(69) 



(70) 



where Q_s = Qgrb/(10 8 Mpc 3 yr 1 ). Since losses limit 
the source distance for fi > 10 19 ' 5 eV to D < fi max ~ 50 
Mpc, one can define an energy Eh such that one expects 
only one source to contribute to the flux at that energy at 
any one time. This corresponds to D b (E b ) = D max and 



E b 



1.6 x 10 21 ZBsQ^h^ 2 



-V 

50 Mpc J 



eV 



(71) 



In terms of D b and fi b , the number of sources contributing 
at any time to the UHECR flux from the sphere bounded 
by Dmax may be written as (Miralda-Escude & Waxman, 
1996): 



N S (E,D< fi max ) 





5 


- E ■ 


D b (E)_ 







(72) 



The distance D b (E) (by definition) corresponds to the dis- 
tance separating the local volume from which the UHECR 
flux at given E will be intrinsically bursty and the extended 
volume beyond where the spread in arrival times is large 
and the fluxes from individual sources overlap, resulting in 
a quasi-steady UHECR flux at Earth. In a similar way, the 
energy E b is the energy above which the UHECR flux should 
become strongly time-dependent, depending on the chance 
occurrence of an UHECR production event within a distance 
Db(E). 

The requirement, implied by observations, that some 
iV s > 10 2 independent sources contribute to the UHECR flux 
above Egzk ~ 10 19 ' 5 eV corresponds to 



D b (10 19 ' 5 eV)< D ' 



: 20 



100e 



D n 



50 Mpc 



Mpc . 



(73) 



This implies a lower limit on B rm s, assuming £ co h ~ 1 Mpc 
and protons (Z = 1): 



>2.0 x 10" 



Q 



-1/2 



D n 



50 Mpc) 



-5/2 



(74) 



In the non-evolving GRB scenario (Q-s ~ 1) this is well 
within the observational constraints on firms. These esti- 
mates correspond to a typical delay time for particles from 
that volume of order 



(tdd) (10 19 - 5 eV,D b ) m 4.9 x 10 3 Qll ^ 



50 Mpc ) 



If the local GRB rate per unit volume is indeed as low as 
proposed by Wijers et al. (1998) in an evolving scenario 
without beaming, Q-s ~ 10~ 2 , the random IGM magnetic 
field would have to be close to the upper limit implied by 
observations: firms — 10 G. In that case one has E b ~ 
3 x 10 20 eV, and the most energetic cosmic ray ever observed 
(£~ 3 x 10 20 eV) is right in the energy range separating 
the quasi-steady flux at lower energies and the intrinsically 
bursty (intermittent) flux at higher energies. 

At energies E > fib, impulsively produced particles ar- 
rive in bursts separated by some 2 x W 2 QZl yr, with the rel- 
ative burst duration scaling as (tdei) /T(_D ma x) = [E/E b ]~ 2 . 
It is therefore unlikely that many particles will be observed 
at any time at energies exceeding ~ 10 21 eV, with E 3> fi b , 
even if they are produced at the source and survive in tran- 
sit. This is especially true in the evolving GRB scenario, 
where fi b ~ 1.6 x 10 20 B- 9 £ eV. Given that some 100 
sources contribute to the flux around 10 19 ' 5 eV, this im- 
plies that the flux should become intrinsically bursty at an 
energy E b ~ 10 20 ' 5 eV. The observation of a single particle 
at 3 x 10 20 eV is therefore within the possibilities of this 
scenario. 

Present experiments have too small a collecting area to 
experimentally check the prediction of a quasi-steady flux 
below fi b and an intrinsically bursty flux above fi b with a 
'duty cycle' decaying as fi~ 2 . The proposed Auger observa- 
tory, which improves the available collection area by more 
than an order of magnitude, should be able to settle this 
question. 

The effect of the delays on the spectrum from a bursting 
source is shown in Figures 12, 13 and 14. Figure 12 shows the 
accumulated arrival spectrum due to a monochromatic burst 
of particles at fi ln j = 10 21 eV for different integration times 
after the arrival of photons from the same source. Particles 
at lower energies accumulate a longer delay, leading to a low- 
energy decline in the spectrum for finite integration time. 

Figure 13 shows a similar calculation, but now for a 
source at a distance of 25 Mpc injecting a power-law spec- 
trum N(E) oc fi~ 2 ' 5 . The assumed IGM field has an ampli- 
tude fi rms = 10" 10 G. Shown are the accumulated particle 
distributions for integration times of 1, 10, 25 and 50 yr 
starting immediately after arrival of photons from the same 
source, and the total arrival spectrum for infinite integration 
time. 

Figure 14 shows the spectrum of particles received in 
five 20-year intervals from a source 25 Mpc distant, with a 
power-law injection spectrum N(E) oc fi~ 2,5 . These parti- 
cles propagate through a random field with a field strength 
fi r = 10~ G with a coherence length of £ cob — 1 Mpc. 
Particles of lower energy are seen with a larger delay, with 
a fairly narrow spread around the mean arrival energy. The 
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accumulated energy distribution 



scaled energy distribution 



t del ^ 1 month, year, decade and 
E ir: = 10 21 eV 
B„ 



rms 

D = 10 Mpe 



10 " 10 G > 'ooh = 1 M P C 




Figure 12. Accumulated particles from a single bursting source 
at 10 Mpc, as a function of the integration time after photon ar- 
rival. The source produces particles at an energy _E; n j = 10 21 eV. 
The assumed intergalactic field has an amplitude B rms = 10~ 
G. The time < de i(max) ranges from one month to one decade. The 
thick line is the total arrival distribution E N(E) for 
For finite integration time there is a sharp decline towards low 
energies due to the strong energy dependence of the delay time, 
(tdel) °c E~ 2 . At the injection energy, the expected delay is 
*dcl — 0.3 year. The finite spread in energy of the particles is 
due to Poissonian statistics of the photon encounter rate in the 
photo-pion production process responsible for the energy losses 
from 10 21 eV. 



spectrum observed at Earth consists of many such contribu- 
tions from different sources with different delays. 

We have simulated arrival spectra that result from a 
distribution of identical bursting sources which produce par- 
ticles with an energy distribution N(E) dE = kE~ s AE. Re- 
sults are shown for s — 2.5. Within a distance D < D[ n = 
100 Mpc, the source distances are drawn from a distribu- 
tion satisfying the scaling law (^) . In addition we calculate 
the quasi-steady background from sources with a distance 
D > _Di n up to a distance Dmax = 2.5 Gpc. This distance 
is so large that, for the assumed field strength in the in- 
tergalactic medium, no particles reach the observer within a 
Hubble time from larger distances (see Eqn. 34) . The sources 
at distances D < Di n are given an 'age' with respect to 
the arrival of photons, distributed uniformly in the range 
< t < (tdel) (E,D), when UHECR detection starts. Dif- 
ferent assumptions about source age have been tried (e.g. 



^ 1 yr, 10 yr, 35 yr, 50 yr and 
N inj (E) = kE- 

B ™ s = 10_10 G, U = 1 M P C 
D = 25 Mpc 




Figure 13. The arrival energy distribution accumulated from 
a source at 25 Mpc after propagation through a 10~ 10 G ran- 
dom field with correlation length ^ con = 1 Mpc. Particles are 
injected with a power-law, N(E) = kE~ 2 ' 5 . The arrival distribu- 
tion is multiplied by a factor <x E 2 T ^ so that an unmodified source 
spectrum corresponds to a horizontal line. The integration time 
ranges from 1 year to 50 year after arrival of photons from the 
same bursting source. The thick line is the total arrival spectrum 
after all particles from the source have been received. 



0.5 (tdel) < t < 1.5 (tdel)) but do not lead to significantly 
different results. The representative energies for determin- 
ing source distribution and age are chosen logarithmically 
equidistant, with E n = E b /(V2) {n ' 1} and n = 13. b igures 
15 and 16 show two representative examples of the resultant 
total spectrum for an assumed field strength B TrnB = 10~ 10 
G. They correspond to two different statistical realizations 
of the the source distribution for D < Dj n . We assumed a 
GRB rate per unit volume of Q_g = 1. At energies below 
E, ~ 10 17 eV, one sees the effect of the transition to spatial 
diffusion as a decline in the spectrum towards lower energies, 
as discussed in Section 2.3. 

These two examples show that there is a significant 'tail' 
of particles above the GZK cut-off energy. However, the flux 
in this tail depends sensitively on the observation window. 
Only a few sources contribute for energies E < Eh, with 
the distance distribution oc D strongly favoring the largest 
possible distances around D — D max . As a result, like in 
the case of continuous sources, the contribution at the high- 
est energies will be dominated by those closest few sources 
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scaled energy distribution 



: t del 0-20, 40-60, 80-100, 120-140, 160-180 yr 
; and for t , , S«= 



" N (E) = /cE- 25 




Figure 14. Particles accumulated in five twenty-year intervals, 
equally spaced in the range < (idol) < 180 y r i from a sin- 
gle bursting source at a distance of 25 Mpc, after propagation 
through a random field with field strength of 10 -10 G and a co- 
herence length £ coh = 1 Mpc. The injection spectrum is a power- 
law, N(E)dE = KE~ 2 - 5 dE. The dominant power-law behaviour 
has been scaled away so that the unmodified source spectrum is 
a horizontal line. The thickest line is the total arrival spectrum 
which would be accumulated over a large time interval after the 
arrival of photons from the same source. 

which in addition have burst at the right time for the parti- 
cles to reach the observer at the present epoch. 



6 DISCUSSION 

We have considered the propagation of ultra-high-energy 
cosmic rays in a disordered but space-filling intergalactic 
magnetic field. We have derived theoretical results for the 
accumulated deflection angle and time delay as a function 
of source distance, and compared these with numerical sim- 
ulations of the random magnetic deflection of UHECRs. 
Depending on the typical magnetic field strength, the de- 
lay for particles originating from a distance of 50 Mpc at 
10 20 eV could range from 10 years for B — 0.01 nG to 
10 5 years if B = 1 nG. We found very good agreement 
between theory and simulations. We have also shown that 
spatial diffusion due to random magnetic cells does not oc- 
cur for distances up to 100 Mpc except for particles with 
energies below ~ 10 19 B-9 eV. Spatial diffusion does limit 



mean flux bursting sources 




Figure 15. The time-avcraged fluxxi? 2 ' 5 due to bursting 
sources, ft is assumed that all sources inject an identical spec- 
trum N(E) AE oc E~ 2 - 5 dE for 10 17 < E < 10 22 eV. The thick 
line in the range 10 17 < i? a rr < 10 19 ' 5 cV corresponds to the 
steady background due to sources more distant than iOO Mpc, 
the burst distance at an energy E\, ~ 10 21 eV for the assumed 
field strength, B r ms = 10~ 1C) G. The decline in this contribution 
below 10 18 eV is due to the fact that particles undergo spatial dif- 
fusion at these energies, and only particles from relatively nearby 
sources can reach the observer within a Hubble time. The sharp 
cut-off at 10 19 ' 5 eV in the flux from distant sources is the GZK 
cut-off due to photo-pion production on the CMWB photons. The 
high-energy part of the spectrum, which peaks around the GZK 
cut-off, is due to bursting sources within 100 Mpc. They are dis- 
tributed according to the N(< D) oc D 5 relation. The mean flux is 
shown for three integration times: t; nt = 1 (thin line), 10 (thicker 
line) and 100 year (thickest line). This flux varies considerably 
above 10 20 eV. 



the background flux from distant sources below the Greisen- 
Zatsepin-Kuz'min cut-off energy so that the UHECR flux 
should decay below an energy of ~ 10 17 — 10 18 eV as the 
sampling volume decreases with energy. 

We have calculated the arrival spectra which result from 
the interaction of UHECRs with the photons of the Cosmic 
Microwave Background Radiation, and demonstrated that 
a significant flux of particles remains above the Greisen- 
Zatsepin-Kuz'min cut-off energy E c ~ 10 19 ' 5 eV for source 
distances up to ~ 50 Mpc. 

We discussed the implications of impulsive production 
of UHECRs in gamma-ray bursts, as proposed by Waxman 
(1995a) and Vietri (1995). We have produced composite 
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mean flux bursting sources 




Figure 16. As in figure 15, but with a different Monte Carlo 
realization of the source distribution N(< D) oc D 5 for D < 100 
Mpc. 

spectra that demonstrate that, just as in the case of contin- 
uous production in powerful radio galaxies, it is likely that 
the highest energy cosmic rays originate from only a few 
sources, located at D < -D max — 50 Mpc. These simulations 
imply that the fact that the GZK cut-off around 3 x 10 19 eV 
is not very pronounced in the observations is probably the 
result of a few UHECR production events relatively close by. 
The Poisson statistics of photo-pion production lead to a sig- 
nificant high-energy tail in the particle flux above the GZK 
cut-off provided enough sources are located within 50 Mpc. 
The observational fact that UHECRs above the GZK cut- 
off energy arrive from directions distributed over the whole 
sky implies, in our view, that the impulsive production sce- 
nario is more likely, given the low space density of powerful 
radio galaxies, which seem to be the only other likely pro- 
duction site of UHECRs except possibly shocks associated 
with large-scale structure formation in the universe, which 
are equally sparse. This conclusion of course only holds if 
UHECRs are hadrons and not some exotic particle that is 
not subject to a significant energy loss mechanism other than 
redshift as a result of the universal expansion. 

In the scenario of impulsive production in gamma-ray 
bursts, the fact that ~ 100 events are seen above 10 19 ' 5 
eV already puts a lower limit on the IGM magnetic field 
strength of B T ~ 0.1 — 1 nG, depending on the GRB rate 
per unit volume. For such a field strength, there should not 
be a significant flux of cosmic rays arriving from outside 



our Galaxy below 10 17 eV, as slow diffusion does not al- 
low such particles to reach us from distant sources within a 
Hubble time. Also, the typical energy scale separating the 
quasi-steady flux at lower energy, due the spread in arrival 
times, and the intrinsically bursty flux at higher energies is 
close to that of the highest energy cosmic ray ever observed: 
E h ~ 10 20 - 5 - 10 215 eV. Our simulations show that the 
resulting flux at Earth in this scenario is already strongly 
variable in time at an energy E ~ 0.2 E^. Therefore any 
interpretation of the observed arrival spectrum in terms of 
a power-law fit, appropriate for continuous sources, is in 
our view hazardous at high energies, as is a firm conclusion 
about the presence (or absence) of a signature of the GZK 
cut-off in the observed spectrum. 

We have made a significant simplification in that we 
assumed that the IGM magnetic field, although randomly 
directed with a coherence length in the Mpc range, fills in- 
tergalactic space homogeneously. A similar calculation could 
be done assuming the magnetic field is concentrated mainly 
in clusters. In that case, one should consider the possibil- 
ity that the cluster fields are so strong (B cl ~ 10~ 7 - 10~ 6 
G, with a possible ordered component), that true spatial 
diffusion of UHECRs could occur in a cluster for particles 
up to an energy of 10 20 — 10 21 eV. This would lead to a 
cluster-scale leaky box model, where UHECRs are confined 
to clusters for a significant time (Rachen, private commu- 
nication), analogous to the situation which prevails within 
the Galaxy for the lower energy cosmic rays. Some calcula- 
tions along these lines have been performed recently by Sigl, 
Lemoine and Biermann (1998) for the local supercluster. In 
that case, most of the UHECRs observed at Earth would be 
from the local group of galaxies, making the question of the 
production site even more pressing. 



ACKNOWLEDGEMENTS 

This research was supported by the NWO-ASTRON foun- 
dation, grant nr. 781-71-050, (AA, YG), the Research Cen- 
ter for Theoretical Astrophysics, University of Sydney, Aus- 
tralia (AA) and the European Union through EEC-HCM 
Contact Nr. ERBCHRX-CT94-0604 and EU-TMR Contract 
Nr. ERBFMRX-CT98-0168. 

REFERENCES 

Achterberg, A. & Kriills, W.M. 1992, A&A, 265, L13. 
Aharonian, F.A. & Cronin, J.W. 1996, Phys. Rev. D 50, 1892. 
Bcrczinskii, V.S., Grigor'eva, S.I. & Zatsepin, G.T. 1975, Astroph. 

Space Sc. 36, 3. 
Berezinskii, V.S. & Grigor'eva, S.I. 1988, A&A, 19, 91. 
Berczinskii, V.S., Bulanov, S.V., Dogiel, V.A., Ginzburg, V.L., 

& Ptuskin, V.S. 1990, Astrophysics of Cosmic Rays, North 

Holland Publ. Co., The Netherlands. 
Biermann, P.L. 1995, Nucl. Phys. B Proc. Suppl. 43, 221 
Bird, D.J. ct al. 1994, ApJ, 424, 491. 
Bird, D.J. ct al. 1995, ApJ, 441, 144. 
Blumenthal, G.R. 1970, Phys. Rev. D., 1, 1596. 
Chodorowski, M., Zdziarski, A. & Sikora, M. 1992, ApJ, 400, 181. 
Davidson, R.C. 1972, Methods in Nonlinear Plasma Theory, Ch. 

8, Academic Press, New York, 
de Lapparent, V., Geller, M.J. & Huchra, J.P. 1986, ApJ, 202, 

LI. 



© 0000 RAS, MNRAS 000, 000-000 



18 A. Achterberg et al. 



Djorgovski, S.G. et al. 1997, Nature, 387, 876. 

Elbert, J.W. & Sommers, P. 1995, ApJ, 441, 151. 

Farrar, G.R. & Biermann, P.L. 1998, Phys. Rev Lett. 81, 3579. 

Gallant, Y.A. & Achterberg, A. 1999, MNRAS, 305, L6. 

Gardiner, C.W. 1983, Handbook of Stochastic Methods, Springer 

Verlag, Berlin. 
Greisen, K. 1966, Phys. Rev. Lett., 16, 748. 
Hayashida, N. et al. 1996, Phys. Rev. Lett., 77, 1000. 
Hill, C.T. & Schramm, D.N. 1985: Phys. Rev. D, 31, 564. 
Kang, H., Ryu, D. & Jones, T.W. 1996, ApJ, 456, 422. 
Kang, H., Rachen, J. P. & Biermann, P.L. 1997, MNRAS, 286, 

257. 

Kewley, L.J., Clay, R.W. & Dawson, B.R. 1996, Astrop. Phys, 5, 
69. 

Kulsrud, R.M., Cen, R., Ostriker, J. P. & Ryu, D., 1997, ApJ, 480, 
481. 

Mannheim, K. & Biermann, P.L. 1989, A&A, 221, 211. 
Miralda-Escude, J. & Waxman, E. 1996, ApJ, 462, L59 
Kronberg, P.P. 1996, Space. Sc. Rev., 75, 387. 
Lesieur, M. 1987, Turbulence in Fluids, Kluwer Academic Publ., 
Dordrecht. 

Metzger, M.R. et al. 1997, Nature, 387, 878. 
Milgrom, M. & Usov, V. 1995, ApJ, 449, L37. 
Miralda-Escude, J. & Waxman, J. 1996, ApJ, 462, L59. 
Narayan, R., Paczyhski, B. & Piran, T., 1992, ApJ 395, L83. 
Norman, C.A., Melrose, D.B. & Achterberg, A. 1995, ApJ, 454, 
60. 

0ksendal, B. 1991, Stochastic Differential Equations, third ed., 

Springer Verlag, Berlin. 
Paczyhski, B. 1986, ApJ, 308, L43. 
Paczyhski, B, 1997, ApJ, submitted 
Rachen, J. P., & Biermann, P.L. 1993, A&A, 272, 161. 
Sigl, G., Lemoine, M. & Biermann, P.L. 1998: MPIFR Preprint 

nr. 774. 

Sigl, G., Schramm, B.N. & Bhattacharjee, P. 1994, Astrop. Phys., 
2, 401. 

Stanev, T., Biermann, P.L., Lloyd-Evans, J. &: Watson, A. A. 

1995, Phys. Rev. Lett., 75, 3056. 
Stepney, S. & Guilbert, P.W. 1983, MNRAS, 204, 1269. 
Sturrock, P.A. 1994, Plasma Physics, Cambridge Univ. Press, 

Cambridge, UK, Ch. 18. 
Subramanian, K. & Barrow, J.D. 1998, Phys. Rev. D, 58, 083502- 

1. 

Takeda, M. et al. 1998, Phys. Rev. Lett. 81, 1163. 

Taylor, G.I. 1921, Proc. London Math. Soc. 20, 196. 

van Paradijs, J. et al. 1997, Nature, 386, 686. 

Vietri, M. 1995, ApJ, 453, 883. 

Waxman, E. 1995a, Phys. Rev. Lett., 75, 386. 

Waxman, E. 1995b, ApJ, 452, LI. 

Waxman, E. & Coppi, P. 1996, ApJ, 464, L75. 

Waxman, E. & Miralda-Escude, J. 1996, ApJ, 472, L89. 

Wijers, R.A.J.M., Bloom, J.S., Bagla, J.S. & Natarajan, P 1998, 

MNRAS, 294, L13 
Yoshida, S. et al. 1995, Astropart. Phys., 3, 105. 
Yoshida, S. & Teshima, M. 1993, Prog. Theor. Phys., 89, 833. 
Zatsepin, G.T. & Kuz'min, V.A. 1966, JETP Lett., 4, 78. 
Zweibel, E.G. & Heiles, C. 1997, Nature, 385, 131. 



APPENDIX A: SCATTERING ON A 
SPECTRUM OF RANDOM FIELDS 

We consider the general case of a broad-band spectrum of 
random turbulent magnetic fields, given by a Fourier expan- 
sion of the form 



B(x) 



d k B(k) e ikx 
(2^)3 B[k) 6 



(Al) 



with (B) = 0. For homogeneous isotropic turbulence the 
Fourier components satisfy the relation 



(B,{k)Bj{k')) = Ml P tJ (2^) 6 8\k + k' 



(A2) 



Here we have employed the standard random-phase approx- 
imation, and defined the projection tensor Py = 8ij — K,iKj 
in terms of the unit vector k = k/\k\ which appears because 
the magnetic field is solenoidal: 



V • B = 



k ■ B(k) = 



(A3) 



The power spectrum B(k) of the magnetic field has been 
normalised in such a way that 



B r 



(\B\ 2 )- / 
Jo 



dkB(k) 



(A4) 



with k = \k\. 

The equation of motion for an ultra-relativistic par- 
ticle yields, by the Kubo- Taylor formula, the following ex- 
pression for the diffusion coefficient describing the propaga- 
tion direction change: 



Ze\ 2 



ds I dk I B(k) 



„ik ■ ns 



(A5) 



We have employed the quasi-linear approximation where the 
unperturbed motion, x(s) = sn, is used to evaluate the phase 
factor in t he integral. To obtain this expression we have used 
Eqn. (A2), defined the solid angle dQk in Fourier space so 
that d^fc = k 2 dfifc dk, and introduced the tensor 



(n x n)i(n x n) 



(A6) 



Choosing n = e z and defining position angles in Fourier 
space by 



k = (sin 8 cos <f> , sin 8 sin <f> , cos 8) 



(A7) 



so that k ■ n = cos 8 = (i, one can perform the integration 
over df2k = d(j> d/i and s. The integration over <f> is simple, 
involving the integral 



> Qij = 2-7T 



1 + A*" 



(5ij - mnj) . 



(A8) 



The integration over path length s only involves the complex 
phase factor: 



ds e lfcMS = ty Stkfj,) — — 5([i) 
k 



(A9) 



The final integration over /i is now trivial, and the resulting 
expression can be written in the same form as for a collection 
of randomly oriented magnetic cells, 



T>ij =T> (8ij —rurij) , (A10) 
where the scalar diffusion coefficient T>o is given by 

" dk B ^ 



m 7T (Ze\ 2 r 



k 



(All) 



The low wavenumber cut-off fc m i n ~ 2n/r s (E) must be intro- 
duced since cells with a size exceeding the particle gyration 
radius will lead to a large-angle deflection in a single cell. 



Consider a power spectrum of the form 



B(k) = B 2 4 



(k£ c 



1 + (k£ c ) 



a+f3 



(A12) 
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with a,/3 > 0. We assume that r s (E) 3> i c for simplicity. 
This spectrum grows as B(k) oc k a for long wavelengths 
(k£ c <§C 1), and decays as B(k) oc k~^ at small scales 
{klc 3> 1). Standard Kolomogorov-type MHD-turbulence, 
where most energy has been fed into the turbulence at a 
scale £ c , predicts a ~ 4 a nd /3 ~ 5/3 (e.g. Lesieur, 1987). 
Substituting this in Eqn. (All), one finds 



V 



n (ZeB \ 2 

8 Cap \—) 4 ' 



where the constant C a p is given by 



a + f] 



a + P 



(A13) 



(A14) 



This result is completely analogous with the corresponding 
expression in the simple case of randomly oriented magnetic 
cells with identical field strength Bo, defining the effective 
coherence length as £ co h = (37rC a ^/4)£ c . For the canonical 
numbers a — 4 and /3 = 5/3 one finds C a p — 0.695 and 
4oh ^ 1.644. 



APPENDIX B: 
DELAYS 



DEVIATION ANGLES AND 



The scattering of UHECRs on randomly oriented magnetic 
cells can be described as a Stochastic Differential Equation 
(SDE) of the Ito form (e.g. 0ksendal, 1991; Gardiner, 1983). 
Employing the summation convention in what follows, the 
stochastic equation of motion for the three components of 
the unit vector n = (ni ,712,^3) along the direction of flight 
reads: 

dm = -2© m ds + P i3 (n) dWj . (Bl) 

Here Pij(h) = 5y — mrij is the projection tensor onto the 
plane perpendicular to n. This tensor satisfies P — P^ and 
PijPjk Pik- 

The quantity dW = (dWi , dW 2 , dW 3 ) is a three- 
component Wiener process satisfying: 



(dWt) = , (dWi dWj) = ds & k 



(B2) 



Here the average (• • •) denotes an ensemble average which 
commutes with integration and differentiation with respect 
to the path length s. The first term in Eqn. ( [Bl[ ) corresponds 
to 'dynamical friction' due to the n-dependence of the diffu- 
sion tensor, and the second (stochastic) term describes true 
scattering. Note that the tensor P projects dW onto the 
plane perpendicular to n, a property employed in the nu- 
merical realization of this process. In the Ito interpretation 
of SDEs the coefficients in a SDE are non-antictpating, i.e. 
they are independent of the stochastic increments dn;. 

From these definitions one can immediately derive the 
following relation for the stochastic integrals n 4 (s): 



d {rij) 
ds 



-2O {ni) 



(ni) (s) = {m) e 



-21> s 



(B3) 



We use a subscript to denote initial values at s = 0. The 
Ito rules for two stochastic integrals X and Y (e.g. 0ksendal, 
1991), 



d(X ■ Y) = dX • Y" + X • dY + dX ■ dY 



(B4) 



together with the relations 
(ds) 2 = 

dW t ds = ds dW z = 

dW t dW 3 = ds Sij , 
lead to the following equations: 
d(ni?ij) = 2T>o (Sij — 3 UiTij) ds + 



(B5) 



(B6) 



+ V2T% {n t P jk dW k + ■ 
Taking an ensemble average one finds an equation for {nitij} 
d{mnj) 



ds 

Integration yields 

(n t nj) = (ni%> < 



= 2Vo Sij — 6X>o {niiij) 



-6T> s 



+ -7T 1 



-6X> s 



(B7) 



(B8) 



It is easily checked that these solutions represent diffu- 
sion of the direction of flight with diffusion coefficient T>ij = 
T>o Pij(n). Making a first-order expansion for T>qs 1, as- 
suming that all members of the ensemble start with the 
same (statistically sharp) initial conditions mo, one finds 
for Am = m — (m)- 



{Am An,) 
2s 



T>o (8^ — morijo) = T>ij(no) 



(B9) 



Next we consider the the distance integrals along the coor- 
dinate axes, 



Xi(s) 



ds' m(s') 



(BIO) 



Taking the average, using equation (B3), one has 



(Bll) 

From the rule (|B4|), together with dxi = m ds, it follows 
that 

d (xiUj) = [mrij ~ 2£> XiUj) ds + V2A) XiP jk dW k . (B12) 
An average yields the differential equation for {xiUj}: 

d (XiTlj) 



ds 



(mrij) - 2T> (xiUj) 



(B13) 



Integrating this equation using Eqn 

/ \ 1 \ { n i n j)o ~ 3 -2T> s /, 

<aWi>W = e (1 



one finds: 

-4-D s\ 



(B14) 



+ 



6V 



-2I>o s 



Note that this expression is symmetric in i and j. In a similar 
fashion one can derive an equation for (xiXj): 

d(xiXj) 



ds 



2 (xiUj) 



(B15) 



where we used the above symmetry. Direct integration 
yields: 
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+ 



3©„ 



-2"D s 



H 1 



(B16) 



1 

2% 



These expressions can be used to calculate time delays 
and deflection angles in the small-angle limit where T>os <iC 
1. Assuming that all particles start from the origin along 
the z-axis so that n(0) = (0 , 0,1), one finds the following 
expressions to first order in T>q: 

(n x ) = (n y ) = , (n z ) = 1 - 2T> s ; 



(nl) = (n%) = 2V s , (n 2 z ) 



(x) = {y)=0 , (z) = s - V s 2 ■ 



(x 2 ) = (y 2 ) = IV s 3 , (z 2 ) = s 2 ~2V s 3 ; 



(B17) 



(xn x ) = (yn y ) = V s 



(zn z ) = s — 3T>os 



The distance r of the particle from the o rigin follows, for 
|a;|,|3/| <§C \z\ — s, from expanding r = \J x 2 + y 2 + z 2 . To 
lowest order in V s 1, the resulting average follows from 
the above expressions: 



(r) M 



(- 2 ) + <S/ 2 ) 



2s 



(B18) 



It, 2 



This determines the average delay time of charged UHECRs 
with respect to ballistically propagating particles (photons) 
originating from the same source at a distance r, to first 
order in Vqs: 



c (tdei) ^ | T> r 2 



(B19) 



The unit vector f — (x , y , z)/r along the line of sight 
from the source to the particle position r(s) at distance s 
makes an angle 6 = cos - 1 (z/r) with the initial flight direc- 
tion. For symmetry reasons, assuming isotropic production 
of particles at the source, one has (6) = 0. The rms value 
#rms = \J (0 2 ) follows from expanding z/r ~ 1 — (x 2 +y 2 )/2z 2 
together with cos 9 ~ 1 



lp2. 
2 r 



{?) = 



(* 2 ) + (y 2 ) 



= § Do » . 



(B20) 



In a similar fashion, one can calculate the rms angle between 
the line of sight from the source to the observer and the flight 
direction from 



n x x + n y y + n z z 



(B21) 



Expanding the cosine and using the above averages, one 
finds: 



(a' 2 ) 



V s . 



(B22) 



This determines the rms angular spread in arrival directions 
around the direction to the source as measured by an ob- 
server at a distance r w s — ct from the source. 



One can also calculate the long-time behaviour for 
T>qs 3> 1. In particular one finds that almost all memory 
of the initial conditions is lost: 



(n.) 

(x) = {y) = 

8ij_ 
6£> 



(n?) 



i . 

3 ' 



(B23) 



(XiUj) 



(XiXj) 



3X»o 



The last relation corresponds to isotropic spatial diffusion, 
where the distance r to the source increases as 



s 



2 



(B24) 
(B25) 



where the spatial diffusion coefficient is 

Here we use s — ct for ultra-relativistic particles. The 
corresponding mean-free-path, defined by K, = iA m f p c, is 
A m fp = 3/2Do- 
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